/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms to convert to the FK5 standard reference frame
Created: PJN / 29-12-2003

Copyright (c) 2003 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/
using System;

namespace AAPlus
{
	public static class FK5
	{
		public static double CorrectionInLongitude(double Longitude, double Latitude, double JD)
		{
			double T = (JD - 2451545) / 36525;
			double Ldash = Longitude - 1.397*T - 0.00031*T*T;

			//Convert to radians
			Ldash = CoordinateTransformation.DegreesToRadians(Ldash);
			Longitude = CoordinateTransformation.DegreesToRadians(Longitude);
			Latitude = CoordinateTransformation.DegreesToRadians(Latitude);

			double val = -0.09033 + 0.03916*(Math.Cos(Ldash) + Math.Sin(Ldash))*Math.Tan(Latitude);
			return CoordinateTransformation.DMSToDegrees(0, 0, val);
		}

		public static double CorrectionInLatitude(double Longitude, double JD)
		{
			double T = (JD - 2451545) / 36525;
			double Ldash = Longitude - 1.397*T - 0.00031*T*T;

			//Convert to radians
			Ldash = CoordinateTransformation.DegreesToRadians(Ldash);
			Longitude = CoordinateTransformation.DegreesToRadians(Longitude);

			double val = 0.03916*(Math.Cos(Ldash) - Math.Sin(Ldash));
			return CoordinateTransformation.DMSToDegrees(0, 0, val);
		}

		public static Coordinate3D ConvertVSOPToFK5J2000(Coordinate3D val)
		{
			Coordinate3D result = new Coordinate3D();
			result.X = val.X + 0.000000440360 * val.Y - 0.000000190919 * val.Z;
			result.Y = -0.000000479966 * val.X + 0.917482137087 * val.Y - 0.397776982902 * val.Z;
			result.Z = 0.397776982902 * val.Y + 0.917482137087 * val.Z;

			return result;
		}

		public static Coordinate3D ConvertVSOPToFK5B1950(Coordinate3D val)
		{
			Coordinate3D result = new Coordinate3D();
			result.X = 0.999925702634 * val.X + 0.012189716217 * val.Y + 0.000011134016 * val.Z;
			result.Y = -0.011179418036 * val.X + 0.917413998946 * val.Y - 0.397777041885 * val.Z;
			result.Z = -0.004859003787 * val.X + 0.397747363646 * val.Y + 0.917482111428 * val.Z;

			return result;
		}

		public static Coordinate3D ConvertVSOPToFK5AnyEquinox(Coordinate3D val, double JDEquinox)
		{
			double t = (JDEquinox - 2451545.0) / 36525;
			double tsquared = t*t;
			double tcubed  = tsquared * t;

			double sigma = 2306.2181*t + 0.30188*tsquared + 0.017988*tcubed;
			sigma = CoordinateTransformation.DegreesToRadians(CoordinateTransformation.DMSToDegrees(0, 0, sigma));

			double zeta = 2306.2181*t + 1.09468*tsquared + 0.018203*tcubed;
			zeta = CoordinateTransformation.DegreesToRadians(CoordinateTransformation.DMSToDegrees(0, 0, zeta));

			double phi = 2004.3109*t - 0.42665*tsquared - 0.041833*tcubed;
			phi = CoordinateTransformation.DegreesToRadians(CoordinateTransformation.DMSToDegrees(0, 0, phi));

			double cossigma = Math.Cos(sigma);
			double coszeta = Math.Cos(zeta);
			double cosphi = Math.Cos(phi);
			double sinsigma = Math.Sin(sigma);
			double sinzeta = Math.Sin(zeta);
			double sinphi = Math.Sin(phi);

			double xx = cossigma * coszeta * cosphi -sinsigma*sinzeta;
			double xy = sinsigma * coszeta + cossigma * sinzeta * cosphi;
			double xz = cossigma * sinphi;
			double yx = -cossigma * sinzeta - sinsigma * coszeta * cosphi;
			double yy = cossigma * coszeta - sinsigma * sinzeta * cosphi;
			double yz = -sinsigma * sinphi;
			double zx = -coszeta * sinphi;
			double zy = -sinzeta * sinphi;
			double zz = cosphi;

			Coordinate3D result = new Coordinate3D();
			result.X = xx * val.X + yx * val.Y + zx * val.Z;
			result.Y = xy * val.X + yy * val.Y + zy * val.Z;
			result.Z = xz * val.X + yz * val.Y + zz * val.Z;

			return result;
		}

	}
}
